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A molecular dynamics (MD) simulation of planar Couette flow is presented for the minimal channel 
in which turbulence structures can be sustained. Evolution over a single breakdown and regeneration 
cycle is compared to computational fluid dynamics (CFD) simulations. Qualitative similar structures 
are observed and turbulent statistics show excellent quantitative agreement. The molecular scale 
law of the wall is presented in which stick-slip molecular wall-fluid interactions replace the no- 
slip conditions. The impact of grid resolution is explored and the observed structures are seen 
to be dependant on averaging time and length scales. The kinetic energy spectra show a range 
of scales are present in the molecular system and that spectral content is dependent on the grid 
resolution employed. The subgrid velocity of the molecules is compared to spatial averaged velocity 
using joint probability density functions. Molecular trajectories, diffusions and Lagrangian statistics 
are presented. The importance of sub-grid scales, relevance of the Kolmogorov lengthscale and 
implications of molecular turbulence are discussed. 


I. INTRODUCTION 

Molecular and turbulent motions are traditionally thought to exist at very different scales, separated by 
three or more orders of magnitude. As a result, turbulence at the molecular scale has not been the focus 
of significant stud 3 ^^, in large part due to the prohibitive molecular system sizes required. However, as will 
be shown in this work, through careful choice of viscosity and geometry a turbulent molecular simulation is 
possible with current computing resources. 

The traditional simulation tool for turbulent flow is continuum computational fluid dynamics (CFD)P1. 
CFD models the flow of fluids away from hydrodynamic equilibrium, using the Navier-Stokes equation. 
The Navier-Stokes equation is based on Newton’s law, with additional assumptions, namely: a Newtonian 
framework in an inertial reference frame, thermodynamic equilibrium (or quasi-equilibrium), fluid isotropy, 
stress tensor symmetry and sometimes incompressibilit}^. Only with the advent of computers have general 
numerical solutions of the Navier-Stokes equations been possible. With these numerical solutions, complex 
and apparently random behavior is the norm, a phenomenon known in fluid mechanics as turbulence. Tur¬ 
bulence is a feature of almost every case of engineering interest. Greater understanding of this phenomenon 
would facilitate improvements in predictive capabilities and advances in engineering design. 

Non-equilibrium molecular dynamics (NEMD)Pis an alternative approach to fluid modeling which requires 
only Newton’s laws and a form of inter-molecular potential. As a more fundamental model, molecular 
dynamics (MD) captures a much greater range of phenomena with no additional assumptions. Examples 
include: arbitrarily strong shock wave^, bubble nucleation, phases change and co-existenc^, moving three- 
phase contact line j^ exact energy conservatioiP, visco-elasticity or memory effect^ and detailed solid- 
liquid interfaces'. The price for this generality is in computational overhead, limiting accessible systems 
to the microscale. Despite this, many complex fluid dynamical phenomena have been reproduced at the 
molecular scale, with excellent agreem ent to the continuum. These include Poiseuille and Couette floytSl^ 
vortex shedding on a cylindeJS^ or plat^SSHH^ Taylor-Couette rollJSSHIll^ the Rayleigh-Taylor instab ilit}^^!!!! 
two dimensional fluid mixin^SHl ^ell as turbulence like vortices due to shock waveP^I. xhe fundamental 
nature of MD means it is uniquely placed to provide insight into the mechanisms of turbulent transport. 
Figure [l] shows a schematic of the spectral energy cascade, believed to be a central feature of turbulent flovP. 
Large scales eddying motions break down into smaller scales until at some minimum length scale the coherent 
motions dissipate to heat as a consequence of viscosity. Energy is conserved in an MD simulation and the 
viscous heating is simply the change of coherent velocity to non-coherent molecular fluctuations. Molecular 
simulation requires no grid and the minimum flow scale is dictated by the underlying configuration of the 
molecules themselves. As a result, molecular modeling of turbulence has the potential to provide a complete 
picture of the energy cascade and link thermodynamic concepts such as pressure and temperature along with 
viscosity to hydrodynamic fluctuations. The range of application of various fluid modeling methodologies 
are included on Figure Reynolds Averaged Navier Stokes (RANS) is the most commonly used industrial 
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FIG. 1: Schematic of the turbulent energy cascade and validity of modeling, including Molecular 
dynamics (MD). Moving left to right, each successive model contains more physical scales of fluid motion 
but is limited by computational cost to simulation of a smaller maximum lengthscale. 


solution to fluid problems and models only the largest scales. Molecular simulation can be viewed as a 
refinement to direct numerical simulation (DNS), where sub-dissipative scales are modelled with no filtering 
of thermodynamic ‘noise’. This is analogous to Large Eddy Simulation (LES) which filters the smaller 
turbulent motions of a DNS, replacing them with a closure model to approximate the details. Erom this 
viewpoint, viscosity and pressure in a continuum solution are simply a closure model to replace the filtered 
molecular detail. There is therefore a clear benefit of MD to industrial fluid mechanics, especially in the 
rapidly developing field of nano and micro fluidics. As with DNS simulation, MD can be used to refine and 
improve the underlying models which are essential to industrial CED simulation. 

In this work, we aim to explore the minimum scale at which turbulence breaks down into viscous heat. We 
reduce the syst em siz e until it contains only the simplest turbulent flow - known as the minimal flow-unit 
in the literatur^^^^. The smallest turbulent flow-unit is found in planar Couette flow, where a viscous 
fluid is driven by two infinite plates sliding in opposite directions. This simulated case is an example of 
wall-bounded turbulent shear flow and, due to the low Reynolds number, only models flow in the near-wall 
region. The flow is not fully turbulent and as a result, the observed flow structures are well organized and a 
repeating cycle of streak ‘busting’ and reformation is observed. In a recent summary, Jimenez^^ concludes 
that available evidence strongly suggests minimal flow units are representative of dynamics in real large scale 
turbulent flow. Experimental evidence also shows that the details of near-wall bursting are identical for a 
range of Reynolds number^^. The minimal channel therefore possess the main attributes of turbulent flow, 
namely that it is non-laminar, three dimensional, stochastic!^ and shows significant and irregular variations 
in both space and tim^. More importantly, due to its minimal size, it is computationally tractable using 
molecular dynamics. 

The first part of the work introduces the modeling methodologies of computational fluid dynamics and 
molecular dynamics. A parameter study is then presented to determine the optimal density and temperature 
to enable a computationally feasible turbulence simulation. Analysis of the flow is performed by identifying 
identical time evolving turbulent structures in both the molecular and continuum systems. The streak 
breakdown and vortex regeneration of the minimum channel flow are demonstrated in the molecular system. 
The phase space evolution is compared for the two systems showing the molecular system is behaving in a 
similar manner to the continuum and that the dynamics are non-laminar. This is then followed by a range of 
standard turbulence statistical analysis. The existence of turbulence at the molecular scale is demonstrated 
through the agreement between the molecular and continuum solutions. Next, a revised law of the wall is 
presented including the molecular stacking region which is unique to molecular simulation. The evolution of 
various forms of energy is plotted and the interchange between hydrodynamic and thermodynamic energies 
is discussed. The impact of different grid-resolutions on the kinetic energy spectrum is explored. Joint 
probability density functions and Lagrangian statistics are employed to gain insight into the sub-grid motions 
of the molecules. This work ends with a discussion of the minimum scales of turbulence and the importance 
of MD as a tool in gain insight into fluid dynamical flows. 
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II. METHODOLOGY 

In this section, the continuum model is presented first to outline how the minimal Couette channel is 
commonly simulated. An initial condition is obtained, which is used in both the continuum and molecular 
study. Next, the molecular dynamics model is presented and a method to obtain the required viscosity 
(and Reynolds number) through selection of temperature and density is discussed. Finally, the setup and 
geometry of the molecular minimal channel case is outlined. 


A. Computational Fluid Dynamics 


Computational fluid dynamics (CFD) is the numerical simulation of fluid motion with dynamics obtained 
from Newton’s law. The equations for an Eulerian control volume V, 


dt 




( 1 ) 


where p is the fluid density, u is the fluid velocity, 11 the pressure tensor and Fbody represents external body 
forces. This can be written at a point in space by applying the divergence theorem, assuming a continuum, 
taking the limit of zero volume and neglected body forces. 


d _ . 

—pu= —V • [puu + nj 


( 2 ) 


By assuming a stress strain closure relation (with stress isotropy, symmetry and thermodynamic equilibrium), 
the pressure tensor, 11, can be approximated in terms of scalar pressure P and velocity, 

U = PI-p [Vu + {Vuf - (V • n) I] - A (V • n) I. (3) 

using Stokes hypothesis A + 2/3/i = 0 and fluid incompressibility, V • n = 0, the Navier-Stokes equations 
are obtained. 


^ + (m-V)m = -VP+TvV (4) 

ot Re 

where the formulas in Eq. 0 are non-dimensionalised by the channel half height h, wall velocity Uw = ±1 
and the kinematic viscosity i'cfd — Mcfd / Pcfd with = hU^ / ^cff> • ChannelflovP^, a spectral solver 

is enmloyed to numerically approximate the solution to the incompressible and isothermal Navier Stokes 
Eq. Q. Channelflow’s numerics are based on Canuto^^, and solve the Navier Stokes equations with a 
Chebyshev-tau algorithm in the primitive variable formulation. To simulate planar Couette flow, the y 
boundary conditions model counter-sliding walls using the no-slip condition, u{Ph) = Pllw. The wall 
boundary conditions and incompressibility are enforced using the tan correctiorP^ with periodic boundaries 
in the x and 2 ; directions. The velocity field is defined on 256 Gauss-Lobatto points in y with 64 by 64 
uniform points in the x and z directions. Chebyshev basis functions are employed in the wall normal 
and Eourier components in the spanwise and streamwise directions. The rotational form of the non-linear 
advection term is used with 2/3 dealiasing. Time advancement is implemented using third order semi- 
implicit backwar ds differencing. A Reynolds number of four hundred is chosen in line with the existing 
literatur ci^^ l ^^ l ^^ l, supported by experiment^ which show turbulence for Reynolds number as low as 360. 
The channel size is chosen to be x Ly x Lz = l.lbnh x 2h x 1.27rh in units normalised by the channel 
half heighfP^. Channelflow is ideally suited for this geometry and has been used extensively to simulate the 
minimal flow unifP^. 

The Reynolds number for the minimum channel is sub-critical and as a result, will not naturally transition 
even in the presences of significant perturbation^^. As in the original work of Hamilton et al.“ turbulent 
flow is obtained by initializing a high Reynolds number case with random perturbations and reducing the 
Reynolds number. This process is used here, with Re = 1000, 700, 500 and finally 400. The CED field at 
Re = 400 was then run for 2000h/I/ time units to ensure any artifacts of the initialization were washed out 
and that turbulence is sustained. The flow field was then saved to provide the initial velocity condition for 
both the MD and CED runs. Channelflow was started from the same initial condition and was run for a 
single flow through time, lOOh/I/, with detailed files collected to allow comparison with the MD solution 
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over the same time. As the continuum minimal channel is driven by the wall shear, the regeneration cycle 
appears to occur indefinitely. A further run over 35 regeneration cycles (3500h/t/ flow though times) was 
then carried out to collect well resolved statistics. 

Reynolds decomposition is employed in order to analyze the evolution of the flow. The mean velocity, 
n, in the channel is defined by averaging over homogeneous directions in both space and time. For some 
arbitrary quantity A, 


A{y) = - :j—— / / / A{x,y,z,t)dtdxdz, (5) 

'^CFD^xJ^z Jo Jo Jo 

where is the entire simulation time. Using the definition of mean velocity from Eq. the velocity 

perturbation can be defined as u' = u — u 


B. Molecular Dynamics 


Molecular dynamics involves the solution of Newton’s law for every molecule in an N-molecular system. 


miVi = Fi, 


(6) 


where rrii is the mass of molecule i, is its acceleration and the force it experiences, F^, is due to the sum 
of all intermolecular interactions, /ij = Here ^ij is the intermolecular potential 

function. 


= 


4e 

0 , 


(A)"-(A) 


( 7 ) 


nj > Vc 


where £ is the molecular length scale, e the energy scale, rrii the mass and = Vi — Vj the inter-molecular 
separation. The pairwise Lennard-Jones potential, Eq. Q, was used with a cutoff of Tc = 2^, which is 
efficient while retaining much of the essential physic^^^^. The Verlet time integration was employed due to 
its excellent energy conservation propertieJ^. The software used in the investigations is designed especially 
for the NEMD styl^roblem presented here, highly optimized for efficient large scale parallel computing and 
thoroug hly ve rifie(P^fH4J^ xhe setup for planar Couette flow is well established in the non-equilibrium MD 
literatur J^ ^ l ^^ l T his includes tethered wall moleculeP^ with temperatures controlled using the Nose Hoover 
t her most afP^l^. The equations of motion for the wall atoms are given by. 


Pi 

- T 

rrii 


Pi =Fi^ - CPi, 

= no + Gkerl) , 




1 

Qi 


' N _ _ 

y^Pn'Pn 


.n=l 




STo , 


(8a) 

(8b) 

(8c) 

(8d) 

(8e) 


where Vi = ri is the molecular velocity; Pijmi is the peculiar velocity in a reference frame moving at the 
wall speed Uw with Ux the unit vector in the x direction; To is the wall thermostat setpoint; = Vi — ro 
the departure of a wall atom from its tethering site ro and is the heat bath mass for the Nose Hoover 
thermostat. The strength of wall tethering is based on the work of Petravic and Harrowell^^ with coefficients 
of /c4 = 5 X 10^ and ko = b x 10^. Only the walls are thermostatted in order to minimize unphysical effects 
in the dynamics of the fiuicP^. 

The Irving and Kirkwood method^, together with the fiuid mechanics concept of a control volume can 
be used to express a molecular system in the same form as the continuum system. The mass in a control 
volume is. 


[ pdV=ly^midi 

Jv \ 


( 9 ) 
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where is defined in terms of a combination of Heaviside functionals so that = 1 when a molecule i 
is inside a control volume and= 0 when i is outsid^^. The angular brackets denote an average. The 
momentum is, 


J pudV = ^ 


( 10 ) 


With Vi=r denoting the laboratory velocity. From the momentum and mass, the so called streaming 
velocity can be defined as. 


u = 



( 11 ) 


The streaming velocity is assumed to be equivalent to the hydrodynamic fluid velocity (continuum Eqs. [TQ. 
By subtracting this from the velocity of a molecule, a definition of kinetic temperature can be obtained. 


T = 




( 12 ) 


where Pilrrii = Vi — u is called the peculiar velocity. Definition of a streaming and peculiar velocity splits 
the laboratory velocity (and associated kinetic energy) of the molecule into a thermal and hydrodynamic 
component. 

The pressure tensor in a molecular system can be calculated by integrating the Irving and Kirkwood^H 
equations over a cubic control volume,^. 


d 

dt 


N 

rriiVidi = 

i=l 


[—puu + n] • dS 


The molecular form of the pressure tensor, 11, is then given by. 



Kinetic Configurational 


(13) 


(14) 


The functional dSi in the kinetic term selects only molecules which have crossed a control volume surface 
during a time period . The kinetic part of the pressure tensor in Eq. (14) is therefore the kinetic theory 
definition; namely a force due to the cumulative effect of molecules bouncing off a container surfac^^. The 
functional dSij selects inter-molecular forces when molecules i and j are on opposite sides of the control 
volume surface (and the inter-molecular force is acting through the surface). The kinetic and configurational 
terms together are the method of planes form of pressure^ localized to the surface of a control volume. 
Mathematically Eq. ( p!4| ) is expressed in a weakened form and is equivalent to the continuum control 
volume equation Q. Note that the term pressure and stress are used interchangeably at the molecular scale 
(where stress is simply negative pressure), although it is natural to speak of the kinetic terms as a pressure 
and the configurational part as a stress. 

The angular brackets should formally denote an ensemble average over many systemJ^, consistent with 
expressing continuum concepts as the expectation of molecular quantified. In practice, provided the average 
phenomena evolves at a slower timescale than the molecular interactions, the angular brackets can be used 
to denote an average over some time interval, , given by. 


{A{x,y,z,t)) = ^— f A{x,y,z,t')dt' (15) 

'^MD J 0 

Despite averaging over the time interval, the left hand side is still a function of a longer timescale. This is a 
distinct difference to the laminar steady state molecular simulation typically studied by NEMD. This results 
in the presence of at least two distinct time scales in a turbulent molecular flows. To obtain the Reynolds 
stress tensor, we start by considering the sum over the outer product of individual molecular velocities. By 
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FIG. 2 : Contour of I/z^^d ^ range of density and temperatures in MD units, interpolated from a 
parameter study of 660 Couette channels (sample points shown by crosses). 


separating out the streaming velocity, we obtain the kinetic part of the molecular pressure tensor in 
(14) and a convective term, puu^ which can be further split into mean flow and turbulent fluctuations. 


Eq. 

i.e., 


'^{rriiriri) = '^{piPi/rrii) + puu + pu'u' 


(16) 


the angular brackets and overbar denote averaging over 1600 and 5.6 x 10^ timesteps respectively. Note that 
for simplicity, decomposition of the scalar pressure, P = P + p', has not been considered here. 

In order to ensure turbulence in the molecular channel, it is desirable to maintain > 400 throughout 

the MD simulation. Determining a molecular viscosity which allows computationally tractable simulations is 
key to the work presented here: the lower the viscosity, the smaller the required MD domain for > 400. 

To this end, 660 independent Couette simulations (domain size 15.8-^ x 15.8^ x 15.8^) were run at different 
density and temperature values with the ratio of stress and strain used to obtain the viscosit}®!. The 
summary of this study is shown as a contour of 1 /z^md Figure]^ A density of 0.3 with low system 
temperatures can be seen to give the lowest viscosity. As a result, a density of = 0.3 and temperature 
T = 0.4 were chosen so 1 /z^md ~ This requires a channel half height oi h = 280.4^ (with 3^ wall, 283.4^ 
total) resulting in around 300 million molecules {N = 291,287,810). The MD domain size was therefore 
1560.4^ X 566.7^ x 1069.9^ and the expected Reynolds number Re = ^ 450. The final choice 

of domain size is somewhat arbitrary, aiming to maximize the simulated Reynolds number, ensure wall 
velocity, density and temperatures are within the range of previous NEMD studieJ^ and still present a 
computationally tractable simulation. 

The wall velocity, = ±1, and consequent strain rate, 7 = 0.0035, are small compared to molecular 
dynamics studied. This is in order to minimize non-linearity in viscosity coefficients, shear heating and 
compressibility effects. 

The temperature T = 0.4 is applied as a thermostat target value to the solid walls which have a density 
= l.O. The initial temperature, T = 0.4, is set by randomly choosing velocities for molecules in the 
face centered cubic lattice. In the liquid region, the molecules are initialized at the same density as the walls 
and removed to get the density of 0.3. On initialization, the initial structure melts and the average domain 
temperature at the start of the simulation drops to T ~ 0.377. As in the continuum solver, a Reynolds 
number of around 400 is sub-critical and a careful choice of initial condition to ensure turbulent flow is 
requirecP^. Having setup the domain, the molecular velocities were adjusted to impose the starting velocity 
field. The same initial velocity field as used in the CFD run, obtained from an initial random perturbation at 
higher Reynolds number and steadily reduced using Channelflow, is applied on a cell by cell basis (uniform 
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64 X 256 X 64 cells with ^ 286 molecules per cell). An adaptation of the algorithm to adjust mean velocity 
on initialization was employecP^. 

The MD simulation was run with a timestep At = 0.005 for 5.6 x 10^ timesteps, corresponding to the 
lOOh/U units required for a single regeneration cycl^^. Further runs for another 3.7 x 10^ timesteps 
were employed to ensure the observed behavior was repeated. This simulation required approximately 280 
thousand CPU hours, run over either 256 cores (8 core Nehalem CPUs), 360 or 720 cores (Westmere 12 core 
processors). Inter-core communication is facilitated by the message passing interface (MPI]P^. The Imperial 
College London computing cluster, CX2 was used with CPU employed instead of GPUs. Current GPU 
random access memories are typically too small to hold the 300 million molecules required for the current 
system. The latest version of GPU based MD codes, such as HOOMrPl, are beginning to allow distribution 
over multiple GPUs (through MPI style communications). Together with increasing GPU random access 
memories, this may soon allow large runs of the type presented here to be routine on desktop computers. 

Due to the work done by the sliding walls in a Gouette flow simulation, heat is added to the MD system, 
which increases the viscosity and decreases the Reynolds number. From the parameter study, the channel’s 
Reynolds number is approximately proportional to the inverse square root of temperature, Re(T) « 300/V^. 
At the initial temperature of T ^ 0.377, the starting Reynolds number is = 489. The arrow in Figure 

shows the change during a single flowthrough time (5.6 x 10^ timesteps) due to heating, with a final value 
of Re^^ = 424. The minimum Reynolds number for which turbulence is sustained. Re = 360 is shown as 
the solid black line in Figure and the GFD Reynolds number, Re^^^ = 400, is shown by a dotted line. 
The choice of parameters therefore ensured that Re^^ was greater than 400 for the regeneration cycle. Even 
during a second regeneration cycle, the temperature increased to T = 0.52 and the Reynold number is still 
Re ^ 416. 


C. Setup Summary 


The length and viscosity in the GFD simulation are matched to the molecular simulation for simplicity. 
The key parameters of both GFD and MD models are presented in table IIC for comparison 


Quantity 

CFD 

MD 

Lx X Ly X Lz 

I.757rh X 2h X I.27rh 

1560.4^ X 566.71 x 1069.9^ 

Simulation 

Gauss-Lobatto in y 

Nmols: 291,287,810 

Size 

64 X 256 X 64 

Uniform Grid: 84 x 198 x 50 

T 

Isothermal 

Variable 0.377 to 0.500 

At 

Variable, max 8 

0.005 

Rve 

28, 000 

8 

Re 

400 

Variable 489 to 424 


Both GFD and MD Reynolds numbers are based on channel half height, h = 280.4^. The MD simulation 
has walls of thickness 3 on both top and bottom, leaving a fluid region of height 560.7-^. The simulation size 
row contains the grid points for the GFD case and number of molecules (with typical averaging bins) for the 
MD case. The averaging time, is used to define the averaging period in the overbar and angular brackets 
definition of mean flow in the GFD and MD respectively. Both simulations write velocity output files at the 
same frequency (every 8 time units). This is a snapshot for the GFD while the MD field is obtained from 
an average of 64 bins in time, each separated by 25 timesteps to ensure decorrelation. 

In the next section, the results of from these two models are compared. 


III. RESULTS 

In this section, the results from the molecular simulation are compared to the same case simulated using 
computational fluid dynamics (GFD). In order to verify the behavior is turbulent; first, instantaneous struc¬ 
tures are shown to be consistent with the GFD solution and the key features of the turbulent regeneration 
cycle are shown to be correctly captured by molecular dynamics. Next, the dynamic behavior of the system 
is demonstrated by evaluating the input energy vs dissipation phase diagram. The evolution of various forms 
of energy are presented and the impact of increasing temperature is discussed. Time averaged statistics are 












8 



0 20 40 60 80 100 

t (h/U) 


FIG. 3: Contour plots of u velocity on the xz plane at the centreline for MD (colors) and CFD (black 
contours with positive (—) and negative (• • •), separated by 0.1). The times of the six contour plots, a) to 
/), are denoted on the plot (CFD (o) and MD (x)) showing the evolution of whole domain turbulent 
kinetic energy k (magnitude of k is 10^ times display value). 


presented and verified, both against the CFD results and studies from the literature. The law of the wall 
is shown with high resolution velocity averaging near the wall to highlight the expected stick-slip behavior 
between the MD wall and fluid. 

Next, molecular results are presented, with particular emphasis on the unique insights offered by this 
discrete nano-scale modeling technique. The impact of the size of averaging volume is explored, by sub¬ 
dividing the domain into different averaging volumes. The velocity energy spectra in the CFD, laminar 
and turbulent MD systems are compared, with a discussion of the possibility of sub-grid motions, viscous 
dissipation and the minimum scale of turbulence. The probability density functions of velocity in a molecular 
system are evaluated showing the range of fluctuations and the relative subtlety of turbulent behavior. 
Finally, a range of Lagrangian style statistics are presented as an example of the powerful insights offered 
by a molecular model. 


A. Instantaneous structures 

In starting from the initial CFD solution, the expected bursting phenomena is reproduced over a whole 
cycle. The turbulent streaks in the MD system can be clearly observed in Figurej^ as well as their breakdown 
and reformation. The MD velocity fields, n, was obtained by averaging the velocity of molecules on a uniform 
grid of 84 x 198 x 50 cell with 64 samples in time each separated by 25 timesteps to avoid correlation^. 
Figure 1^ shows the contours of velocity at the channel centreline in the MD system. The brea kdow n and 
reformation of the streaks is consistent with well documented behavior in the CFD literatur^^^EZl. The 
overlayed black line contours in Figure a) to c) are the results from Channelflow at equivalent times. 
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(a) Isosurfaces of streamwise vorticity, ujx — —0.005 (b) Isosurfaces of spanwise vorticity, Uy — —0.0035 

(red) and uJx — 0.005 (blue). (red) and ujy = 0.0035 (blue). 

FIG. 4: Vorticity isosurfaces at times, a) to /), from figure 


Qualitatively similar turbulent structure are clearly present in both the CFD and MD simulation. The 
molecules model includes both temperature and density variations. 

The compressibility in the MD models would not be expected to be negligible, with speed of sound 
predicted to be approximately 248m/The wall is sliding at a speed of IbOm/s so the the Mach number 
is therefore about 0.65 

It is therefore surprising that there is such good agreement between an isothermal and incompressible 
continuum solution and the molecular model. It appears that despite significant spatial and temporal 
variation in temperature, the coupling to viscosity is weak and does not appear to effect the hydrodynamics. 
Although there is good agreement for the first part of the regeneration cycle, the solutions diverge as time 
progresses and Figurej^d) to /) are shown for the molecular model only. This divergence is mainly attributed 
to the non-linear nature of both the MD and CFD flows. Exact structures observed in two continuum cases 
started from the same initial condition would also diverge over time. This time evolution of the instantaneous 
turbulent kinetic energy k = is shown below the contours in Figure]^ The reformation of 

the streaks from the point of minimum turbulent kinetic energy. Figure [^d) , is clearly seen in the successive 
MD only contours and the plot of turbulent kinetic energy below Figure]^ 

The cycle shown in Figure is an example the so called burstiM phenomena, which includes the 
break down and reformation of coherent structur es in turbulent ffow^. This has been observed in both 
experimentaP^ and numerical studies of turbulenc j^^ l ^^ i This cycle consists of three parts; first the streaks 
become unstable and break apart (streak break-down), next the stream-wise vortices in the flow become 
stronger (vortex regeneration), finally a new set of streaks form from the streamwise vortices (streak 
formation)!^. The streak breakdown and reformation are clearly observed in Figure]^ 

The formation of vortices at the point of breakdown can be observed by looking at vectors of v and w in 
Figure By comparison with the vortices at the end of the simulation (when the streaks had reformed) 
it is clear that the vortices are much stronger at the point of maximum breakdown. Indeed, a secondary 
vortex appears to have formed, which may be a result of the slightly higher Reynolds number at the time 
of breakdown {Re ~ 430) in the MD system. The vorticity isosurface of Fig. show a large increase in 
streamwise vortices at the point of streak breakdown in the minimal channel cycle. 

Figures and suggest that the key stages of the near wall regeneration, or bursting, cycle are repro¬ 
duced in a molecular simulation. This bursting phenomena is commonly believed to be the fundamental 
mechanisms of turbulent energy productioiP^. For larger syste ms, it has been suggested that turbulence 
cannot be maintained without these near wall bursting structureP^^^. As this fundamental characteristic 
of turbulent flow can be reproduced by a purely molecular model, it could provide a new perspective in the 
simulation of turbulence flow. 

As the initial condition used for the molecular system is based on a turbulence CFD field, it is possible 
that the observed behavior is a result of inertia or memory of that initial field. To further test the persistence 
of this MD regeneration mechanism, the molecular channel at t = 5.6 x 10^ and Re = 424 is run beyond 
the single cycle. A similar breakdown is observed at iteration 7.2 x 10^, followed by vortex regeneration and 
subsequent reformation of the streaks. It seems unlikely that such complex behavior could occur simply 
from an initial condition, especially over two distinct and unique regeneration cycles. 
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z 

FIG. 5: UW velocity vectors and streamlines on the yz plane averaged over all x. Contour shows 
temperature variations relative to the mean of the whole domain at that time. Times are chosen before 
streak breakdown, directly after breakdown and after streak reformation. 



In experimental flows below a critical Reynolds number (about 360) the breakdown of the streaks occur 
prematurely and they do not refornP^. To explore this, two smaller MD simulations are run with the 
same initial condition, temperatures and density but smaller domains containing 2, 048, 308 and 11, 935,488 
molecules respectively. The effective Reynolds numbers are Re = 40 and Re =100 and in both cases, the 
streaks and vortices dissipate and a transition to laminar flow is observed. This behavior is consistent with 
equivalent simulations performed using Channelflow at Reynolds numbers of 40 and 100. The appropriate 
MD domain size for a Reynolds number greater than 400 therefore observes at least two regeneration cycles, 
while in smaller domains, the streaks dissipate and the flow relaminarizes. 

The phase space plots of Figure lend further weight to the assertion that the flow remains turbulent 
throughout the simulation. The dissipation is given by, 

D=^ [ {\Vu\^+ \Vv\^+ \Vw\^)dV (17) 

F Jv 

and the mechanical energy input at the wall is. 



The MD values are averaged over 6400 samples separated by 25 time steps each in order to minimize the 
molecular scale fluctuations in the sample data. In a steady state laminar flow, the I and D values are 
equal, denoted by the dotted line in Figure Both CFD and MD are normalized by the analytical solution 
for the steady laminar Couette flow solution. The numerical results for simulated laminar flow are shown 
at the bottom left of Figure The MD laminar solution is larger than unity due to the inherent thermal 
fluctuations (which has a greater impact on derivatives such as strain rates). 

These phase plots are inspired by chaos theory style representations of dynamical system. The many 
degrees of system are reduced to two central features for wall driven shear, namely the rate mechanical 
energy is input into the flow and the rate at which the mechanical energy is dissipated from coherent velocity 
to heat. For the minimal flow unit, the various stages in the regeneration cycle see greater dissipation at 
some times and more energy input at others. As the dissipation should approximately balance the input on 
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I 

FIG. 6: Plot of energy input vs dissipation normalized by the laminar solution. The CFD trajectory over 
35 regeneration cycles (—) and for the single flow through (— •— ) compared to the MD single flow through 
(— •— ) and beyond (-). The laminar solutions in both cases are shown by (x/x) for CFD and MD 

respectively 


average for a cycle, an orbiting behavior is observed in / vs D phase space. The non-linear hydrodynamic 
behavior is similar in both the CFD and MD model, evidenced by the similar orbiting behavior in Figure 
for both systems. 

Note, however, that the MD dissipation is observed to be lower on average than the CFD. The CFD system 
is in a thermodynamic steady state with a constant viscosity while as the molecular system progresses the 
temperature and viscosity continues to change. The dissipation D and input / only measure the energy 
change in hydrodynamic velocity components. Energy added to hydrodynamic components by the wall, /, 
will predominantly change to thermal energy as measured by the mechanism for dissipation D. The two 
therefore appear to still balance in the MD system. The sliding and thermostatting in the molecular wall 
will exchange energy directly with the thermal components of the MD system, but these do not appear to 
significantly impact the I vs D cycle. 

To further explore the energy change in a molecular system, the evolution of various forms of energy is 
considered next. The energy in the system can be divided as follows. 


fCfiow = \ j p\'^?dv, 

(19a) 

1 JV ^ N 

ICrhem = - uf, 

^ i=l ^ i=l 

1 ^ 

^ — 2 ^ ^ '^i \ '^i 1 ~ ICThem T 1C flow-) 

(19b) 

(19c) 

N N 

(19d) 

i=i j=i 


n = ic + T. 

(19e) 


These include the total system kinetic energy, /C, separated into the the flow mechanical energy 1C flow 
of the streaming velocity u (averaged over 64 uncorrelated samples) and thermal energy ICThem- The 






12 


a) b) c) d) e) f) 

0.90 
^ 0.75 
0.60 

0.942 
0.940 
0.938 


0.004 
^ 1 ^ 0.002 
0.000 

0 20 40 60 80 100 

Time 

FIG. 7: Plots of the energy changes as a function of time (h/U) in the molecular system, labelled with the 

cases from Fig[^ 

i) Total energy, H (—), total system kinetic energy, 1C (--), thermal energy, ICThem (* * *) with potential 
energy, T (— ) and flow mechanical energy, 1C flow (--) on second axis, 

ii) Total kinetic energy, /C/H, and potential energy, T/'H, normalized by total energy shown on different 

axis with same range. 

iii) Rate of change of flow mechanical energy, 1C flow {—) and thermal kinetic energy, ICTherm (--) from 
spline derivatives with actual points shown ( ), dissipation, D/Diaminar {—) and energy input by shear, 

Illlaminar (* * * ) ^^ the SCCOud axis. 



4k).. 



potential energy, T, is obtained from the sum of intermolecular potential interactions. The total energy H 
(Hamiltonian) is the sum of kinetic and configurational contributions. 

Figure]^ shows the increase in the various types of energy, Eqs. (19a 19e), in the channel. The total 
system energy, H monotonically increases as the energy added by shearing is greater than the heat removed 
by the thermostatting in the walls. The total kinetic and thermal energy monotonically increases, while the 
mechanical energy in the flow changes based on the regeneration cycle and appear to stay fairly constant. 
There is an apparent decoupling between the thermodynamics and hydrodynamics energies. Most CFD 
simulations do not observe significant heating. In an MD simulation, the impact of heating is far more 
pronounced due to smaller system sizes and much higher shear rates. The explicit energy conservation in 
the fluid and the choice to apply a thermostat to only the wall molecules results in the continual heating 
observed. The heating has an impact of the effective viscosity in the molecular system which increases as the 
simulation progresses. Despite being away from thermodynamic equilibrium, the impact of this changing 
viscosity is not sufficient to prevent the regeneration cycle. 

Figure [^i shows the evolution of normalised system kinetic energy, /C/H, and potential energy, T/H. By 
normalising using the total energy H, the impact of the increasing system energy is removed and it is clear 
that a decrease in kinetic energy of the molecules results in a corresponding increase in potential energy. 
The interchange of potential and kinetic energy is apparently effected by the breakdown and regeneration 
cycle with corresponding peaks and troughs (c./. 19a in Fig[^). Higher potential energy is indicative of 


clustering of molecules. The hydrodynamic bursting cycle governs the flow of molecules and will result in 
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changes in the density (for example inside the streaks). 

Figure [^ii shows the time evolution of flow energy, 1C flow ^ and thermal energy JCThem- The lines in 
the bottom panel Figure are obtained from the derivative of a spline fit to the data (actual points in 
gray). In addition, the input, I and dissipation D as plotted in the phase plot of Figure]^ are included. The 
dissipation, D, closely matches the change in ICTherm (after the initial part of the simulation), suggesting the 
dissipated mechanical energy is transferred directly to the thermal energy. The time evolution of temperature 
is always positive as heat is continuously added to the system, however the rate is decreasing throughout 
the simulation. In the extended run, not shown, by the end of the second regeneration cycle (approximately 
160h/U) the change in thermal energy appears to be zero on average. 

There is a large decrease in the rate of thermal energy production from around the times c) to d), which 
coincides with the streak breakdown and vortex regeneration. It has been suggested that after streak 
breakdown, t he vor tices are re-energised in a complicated mechanism, requiring the interaction of several 
Fourier mode^^^^. The regeneration of the mean flow takes energy from smaller wavelengths and it is 
possible that some of the re-energization of vortices may be at the expense of what would be considered to 
be sub-grid thermal motion. The division between thermal energy and hydrodynamic energy is arbitrary 
in MD; dependent on the spatial and temporal averaging scales. The regeneration of the vortices could 
therefore be at the expense of energy which contributes to either thermal or hydrodynamic energy. A full 
energy analysis of the MD flow would be required to provide insight into distribution of energy in various 
forms. 


B. Time Averaged Statistics 


In this section, we compare the time averaged statistical properties, both to published data and results 
from CFD simulations using Channelflow. Figure [^compares a range of time averaged statistics as a function 
of wall normal position. Panel [^) shows average velocity profiles for both the CFD and MD cases. There 
is good agreement between the streamwise velocity F, averaged over lOOh/U for the MD system compared 
to a very long simulation for the CFD (3500h/I/). The straight line is the analytical solution for fully- 
developed laminar Couette flow. The flow is driven by the shearing due to the wall, which in a continuum 
model results in a repeating cycle at this Reynolds number. The CFD model was run for over 5500 flow 
through times with no sign of re-laminarization in the continuum case. As molecular dynamics is inherently 
energy conserving, shearing from the molecular walls results in a continual increase in the temperature of 
the molecular fluid. Some heat energy is removed at the wall by a thermostat, but the available surface 
to remove heat is insufficient for a domain of this size. The initial and final profile are displayed in Figure 
[^) for the molecular channel. Despite the decreasing Reynolds number, the profile still suggests the flow 
is turbulent. The shape of this profile is due to the streamwise vortices in the flovl^, which as observed in 
Figureare clearly present throughout. For the smaller molecular channels tested at effective Re = 40 and 
Re = 100, the turbulent streaks decay and the flow quickly relaminarizes, returning to the linear profile. 

Root mean square (RMS) turbulence intensities for each velocity component are shown in Figure]^) for 
both the CFD and MD systems. There is good agreement for the three velocity components, despite the 
stick-slip behavior at the walls, compressibility effects, temperature variations and the varying Reynolds 
number in the MD simulation. In addition, results from the CFD literatur^^ at Re = 400 are compared in 
Figure]^) demonstrating good agreement. 

The shear and turbulent stresses are presented in Figure |^c). Again, good agreement is observed between 
continuum and molecular results. The greater value of turbulent shear stress towards the channel center in 
Figure [^c) is attributed to the higher Reynolds number of the molecular simulation. Molecular dynamics 
does not require a viscosity to be defined, with pressure obtained directly from Eq. (14) as a function 
of the molecular configuration and kinetic motions. The molecular kinetic pressure and configurational 
stress contributions are shown in Figure |^c). These MD pressure results are obtained on the same grid as 
the velocity measurements by recording every single interactions and crossings in a 1600 timesteps period. 
During the simulation, 1500 of these periods (2.4 x 10^ measurements in total) are taken at representative 
intervals. The average kinetic pressure and configurational stress sum to a total shear pressure. This total 
shear pressure can also be calculated in both the MD and CFD cases using {1/Re)du/dy. For the CFD case, 
RCcpd = 400, while for the MD case a Reynolds number of Re^^ = 430 gives the best agreement with the 
molecular pressure. An MD viscosity can be obtained from the average over the changing viscosity, 
using temperature and density to lookup values from the data in Figure The viscosity which gives the 
best fit is approximately 10% lower than predicted from Figure This may indicate that the turbulent flow 
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FIG. 8: Comparison plots of channel statistics for MD and CFD: a) mean (—), first ( — ) and last ( — ) 
velocity profile for MD, mean CFD (o) and laminar analytical (—). b) root mean square turbulent 
velocities '^rms ?'^rms 5 '^rms with MD (—,• • •), CFD literaturd^ (□, o, A) and Channelflow (x,x, x) c) 
shear and turbulent stresses with MD/CFD for {1/Re)du/dy (--/□), u'v' (—/o), MD kinetic pressure ( — ), 
configurational stress ( — ) and total pressure ( — ). MD and CFD stresses in plot c) are normalized by MD 

and CFD viscous wall stresses respectively. 


in the molecular simulation has an impact on the viscosity coefficient. 

Another interesting observation is the relative magnitude of the kinetic and configurational parts of the 
shear pressure. Despite the low density of the fluid, the configurational shear is a comparable magnitude to 
the kinetic shear pressure while the direct configurational pressure contributes only 5% of total pressure. 

Another characteristic of turbulent flow is the law of the wall, obtained by expressing the velocity and 
wall normal position in dimensionless units = y/6r and u'^ = u/u^ where Ur=Twlp and 5^ = vjur. 
In dimensionless wall units, the observed behavior is similar for a range of different Reynolds numbers and 
flow geometry. The definition of wall stress, Tyj is ambiguous in a molecular system as a result of stick slip 
behavior in the near wall fluid as well as deformation and stress in the wall itself. The wall stress, is 
evaluated at a distance of 11 away from the wall, the location at which shear stress equals the continuum 
analytical value in the laminar MD channel of the same dimensions. The stress is obtained using the method 
of planes definition, Eq. (14). 

Figure 9a shows the velocity, against logarithmic wall normal position, both expressed in inner 
units. The expected viscous sub-layer, and buffer regions are labeled. The start of the log law region is 
also indicated but is not expected to be significant at such low Reynolds number^^. Of main interest here 
is the molecular stacking region near the wall which is uniquely captured by molecular simulation. A high 
averaging resolution is employed, with 3168 cell in the y direction to explore the near wall detail. This near¬ 
wall molecular stacking in Figure 9a has been widely observed in MD studieJ^ as well as experiment^. 
The fluid solid inter-molecular interactions result in stick-slip behavior and an effective Navier slip length. 
It has been observed that wall roughness in a turbulent flow impacts the form of the viscous sub-layer and 
log law regiorP^. Given the relative size of the wall features to the characteristic scales of the flow, the 
molecular roughness may impact the flow in nano and micro-scale channels. The law of the wall fit with 
surface roughness uses coefficients n = 0.41 and B = 3.2, perhaps suggests the impact of the molecular scale 
on the mean flow. Higher Reynolds number MD simulations would be required to fully explore this. 

The minimal channel models the evolution of turbulent structures in the near wall region. This is where 
the major ity of the production of turbulent energy occur^. In Figure 9b, the rate of turbulent production, 
V = u'^u'adu^/drp is compared for the CFD and MD simulations. There is good agreement between both 
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(a) Law of the wall in reduced units over 3168 y cells 
shown for sum of cell momentum pujpmean (—), cell 
velocity u (logarithmically spaced o), viscous sub-layer 
analytical solutions (—) and log law region (- -). 



Vlh 

(b) Turbulent production, 7^, in the 
MD (—) and CFD (o) case. 
Right: Average density, p (- -) 
and temperature, T (—) at times 
t = {0,3,14, 26,43, 71,100, im}h/U. 


FIG. 9: Law of the wall and temperature profile in channel 


modeling methodologies with similar profiles and maximum production at the same location. The production 
appears to become negative near the channel centre in both cases, perhaps due to the low Reynolds number or 
insufficient statistics from a single flow through time. On the left of Figure the density and temperature 
in the MD channel are also presented. The density is seen to increase near the wall, which may be due 
to molecular stacking near the higher density solid = 1-0). This may also be a consequence of the 

vortices in the flow, which move molecules near to the wall where they are slowed by the interaction with 
the wall. The temperature is seen to be increasing throughout the simulation, while the wall temperature 
remains at T = 0.4 due to the thermostat. There is a Kapitza like jump near the wall which becomes 
more pronounced as the temperature difference between the fluid and solid becomes greater. The initial 
temperature field is flat as only the velocity field is specified from the CFD initial condition. The small peak 
in the initial temperature plot is a shock wave which appears due to initialization of a system between fixed 
walls. This shock wave bounces between the walls and has dissipated by the next displayed time 3h/U. The 
effect of shear heating results in a continually increasing temperature and as a result, an increasing viscosity 
and decreasing Reynolds number. 

In this section, the MD model has been compared to CFD results, showing good agreement for both 
the key stages of the regeneration cycle and reproduction of turbulent statistics. Having demonstrated the 
molecular model is reproducing the minimal flow unit, in the next section the unique insights presented by 
turbulence in a molecular model are explored. 


C. Spectra 


A grid resolution study is performed in this section, with varying spatial and temporal averaging used 
to explore the flow structures at different scales. Figure 10 shows the impact of averaging on the observed 
velocity field at the center of the channel. The far left hand side shows individual molecules, each coloured 
by their own velocity. Moving right, molecules are colored by the average velocity for the cell they are 
located in. The finest grained example has 800 by 1344 cells, with fewer cells employed in even steps until 
the far right velocity field is shown for the 84 by 50 cell resolution. It is apparent that only by employing 
appropriate spatial averaging can the details of the flow structures be observed. To understand the impact 
of grid resolution on the energy content of the fluid flow, a Fourier transform is taken to obtain a spectrum 
of the various scales of fluid motion. 

In Figure El four levels of grid refinement are considered in the MD channel compared to a single CFD 
case, as outlined in table III C| 
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FIG. 10: Plot of ^ 1.4 million molecules in a central y cells of width 2.86^. Coloring is by individual 
molecular velocity on the far left and moving right shows molecules colored by the average velocity in 
increasingly larger cells. The first insert is colored by cell velocities and the next level shows the velocities 

of individual molecules. 


Case 

'^MD 

Grid resolution 

a) 

0.005 

672 X 198 X 400 

h) 

0.005 

84 X 198 X 50 

c) 

8 

84 X 198 X 50 

d) 

32 

84 X 66 X 50 


These are compared to a single CFD case with a xz spectral resolution of 42 x 42 modes. The Kolmogorov 
five thirds law is shown for reference, although due to the very low Reynolds number, no inertial range would 
be expected. Each turbulent MD spectra is compared to the laminar MD solutions (dotted lines in Figure 


11) which have been averaged in the same manner. These laminar solution are run at the same system size 
but with no turbulent initial condition. The turbulent and laminar MD simulations show comparable base 
levels of broadband ’noise’ in the spectra of Figure pTj while only the turbulent case shows low wavelength 
structures. These low wavelength structures exhibit similar shapes and magnitudes to the spectra obtained 
from the CFD solution. One discrepancy is the lower energy in the CFD x spectra relative to the molecular 
system. This may simply be a consequence of using different numbers of cells in the CFD and MD, the 
divergence of the trajectories or the higher effective Reynolds number in the MD simulation. However, this 
may also be a consequence of the greater energy in an MD system due to the sub grid molecular fluctuations. 
This is supported by the results from case a)^ the finest grid resolution, where the entire spectral energy 
is much higher than observed in coarser grids. The level of thermal motion is at much greater energy and 
only two peaks are observed above this. From Figure pR} it is clear that this trend continues with increased 
grid refinement, until at the level of individual molecular velocities, even the streaks are difficult to identify. 
The energy in the streamwise u component, 1C flow i increases with a finer grid at the expense of the thermal 
energy, ICThem^ ns flow energy is split between thermal and streaming. The higher spectral energy in the 
finest grid of Figure pT] suggests that coarse graining can result in lower energy in important flow structures. 
In addition, the process of averaging results in spatial and temporal resolution being lost as the degrees 
of freedom are reduced. However, without coarse graining, it would be impossible to identify the spectral 
content which may be of central importance to the regeneration cycl^^. For example in the 2 : spectra of 
Figure IE the roll structures give a large peak and a range of successive harmonics are clearly observed 
in the CFD solution. As the level of averaging is increased, spectral content previously obscured below 
the thermodynamic fluctuations becomes apparent in Fig El At least two further harmonics, previously 
hidden by the higher frequency thermal ‘noise’ are uncovered. The new peak are at the same wavelengths as a 
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k k 

FIG. 11: The x (left) and ^ (right) axis spectral energy of the streamwise velocity u at the channel 
centreline, a) A fine grained single snapshot in turbulent ( ) and laminar ( ) flow, b) A single snapshot 

on a comparable grid to the CFD, turbulent ( ) and laminar (- -) flow, c) An average of 64 snapshots, 

turbulent (—) laminar (--) flow, d) 6400 snapshots and three y cell spatial average, turbulent (—) and 
laminar (--) flow. The CFD solution (—), with laminar case at limit of machine precision for all k. 
Kolmogorov law (—) E = Cek~^^^ for inertial subrange. 


harmonic observed in the continuum solution. Further low wavelegth peaks in the CFD spectra suggests that 
further harmonics would also become apparent in the MD spectra if even greater averaging was employed. 
A range of flow length scales are clearly present in the MD minimal channel, as apparent in the energy 
spectra of Fig[TTj while the equivalent size of laminar channel does not display any hierarchy of scales. 

One could question where the Kolmogorov microscale, 77 = (z/^/e)^/^, would be located on the spectrum of 
Figure prj Although the derivation of a Kolmogorov lengthscale assumes a much higher Reynolds number, 
consideration of the minimum eddy length scale at which coherent motion is dissipated to heat is a central 
concept in CFD modeling. Both for ensuring a DNS simulation is well resolved and identifying the minimum 
scale of turbulent eddy. Although the minimum scale does not become apparent in the grid resolution study 
of Figuremolecular dynamics may provide insight into the minimum scale of turbulent eddy. If we assume 
that the Kolmogorov scale is below the observed noise floor and impossible to identify, we consider instead 
the minimum physical size this eddy could have in the molecular paradigm. It is possible to view molecular 
liquids as a series of potential wells, created by the cage like configurational structure of a molecular fluicP^. 
These molecular cages would typically see three dimensions motion of the molecule trapped within. This 
provides an absolute minimum scale for a rotational vortex in a molecular simulation, displayed schematically 
on Figure [2 Although these motions are likely different from the large scale vortical motions observed in 
the minimal channel, the interplay between inertial and the smallest scales of motions is of great interest. 
Molecular dynamics is uniquely placed to provide insight into the minimum scale of eddy. 

This section has demonstrated that the choice of grid has a major impact on the observed behavior of the 
system. In order to explore the velocity without requiring the definition of a grid, the probability density 
function are evaluated in the next section. These distribution are compared to the probability density 
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function of the cell averaged velocity. 


D. Probability Density Functions 


The importance of cell based averaging techniques to the observation of flow structures was demonstrated 
in the previous subsection. By averaging over time and space, information is essentially discarded about the 
range of fluctuations below the grid scale. In this section, the distribution of molecular position and velocity 
are explored using probability density functions. The joint probability density function (PDF) is displayed 
for individual molecular velocities and compared to velocity of the spatially averaged cells. The positional 
distribution of individual molecules is also considered in the form of the Radial distribution function. The 
joint PDF is defined for each bin n of width /S.u and height Au, 


v) = P{[n — l]Au < u < nAu and [n — 1] Au < v < nAu), 


( 20 ) 


with normalisation 


, v)dudv = 1. The PDF for all the x and y velocities of the ^ 1.4 million 

yi) in a single y cell at the channel centreline is shown in Figure [12^ The 


individual molecules, Pu^vi^i^ 
distribution of the molecular velocities is almost perfectly Gaussian for both the x and y velocities. This is 
demonstrated by the Ga ussian fit shown on the x projection at the bottom (similar results can be shown for 
y). Also shown in Figure 12a is the joint PDF of Pu,v{u^ v) where u and v are the time averaged {t^^= 64) 


velocity in the 84 x 50 cells at the channel centreline. A total of 338 u and v samples in time are used to 
give comparable statistics to the molecular case (a total of ^ 1.4 million samples). These coarse grained 
cell velocities are shown as a scatter plot in Figure [T^ with the zoom insert showing a contour plot of the 
PDF. The u projection of the CFD velocity PDF is also shown at the bottom of Figure |12a| As there is 
negligible mean flow at the centreline, the average velocity in this top PDF are the perturbation velocities 
u' and v'. The high speed u streak clearly has an an associated positive v velocity while the low speed 
has an associated negative v. The remaining distribution is distorted indicating the effect of the turbulent 
shear stretching. There is a strong sweep and ejection behavior with a slight preference for ejection seen 
in the cell averaged PDF velocities^. The sweep and ejection are though to be the main contribution to 


turbulent energy production and the cell averaged results obtained in Figure |12a| are consistent with this 
observation. The most striking conclusion, however, is that the probability density function of individual 
molecular velocities would simple obscure all of this detail. 

The Radial Distribution Function (RDF) is shown in Figure |12b | and compared to neutron scattering 
experimental data for liquid Argon. The RDF is calculated fronl^ g[r) = {V/N)[n{r)/A'Kr‘^Ar] where n(r) 
is the number of molecules in a spherical shell of size Ar. To improve statistics, the RDF is calculated 
for every molecule in the system. The available experimental data is at p = 0.708 and T = 0.8352 and 
the RDF for the equivalent state point in MD units is shown. The RDF for the state point used in the 
turbulent study, T = 0.4, p = 0.3, is also included in figure |12b| It can be seen that the Lennard-Jones 
model, even for the short (WCA) cut off range used in this work, closely reproduces the structural properties 
of Argon. This reproduction of the RDF also highlights the greatest strength of molecular dynamics; it is 
unique in being the only modeling technique to explicitly capture the underlying atomic structure. The 
fluid’s dynamics are governed by the interplay of molecular kinetic energy with the continual evolution of 
the current configuration of molecules. As a result, the entire history of the molecular structure has a 
continual impact on the current state. Pressure and viscosity are a consequence of the systems evolution 
and it is for this reason, the continuum style stress closure of Eq. (§ is not required in an MD model. 

To explore the evolution of the impact of the molecules history on the flow, the trajectory of individual 
molecules are investigated in the next section. 


E. Lagrangian Statistics 

Lagrangian statistics have shown great promise as a means to understand the nature of turbulent flow^^. 
This involves experimentally following tracer particles, or in numerical studies, adding imaginary tracer 
particles into the turbulent flow. In molecular dynamics the molecules themselves can be considered as 
tracer particles and Lagrangian statistics obtained from their evolution in time. The distinction in molecular 
dynamics is that molecules are both tracers and the flow itself. A key parameter in Lagrangian statistics is 
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Xi 

(a) The joint probability density function for the 
molecular and averaged cell velocities. The bottom 
plot shows the x projection of the molecular 
distribution (—) Gaussian fit (o) and CFD 
distribution (—). 



(b) Radial distribution function for the current case 
p — 0.3 and T = 0.4 (—) with verification at 
p — 0.708 and T = 0.8352 (x), equivalent to 
experimental results for liquid Argon at T = 85iF, 
/5 = 0.02125A“’^ (—). 


FIG. 12: Probability density and Radial distribution functions 


the structure factor, 


1 ^ 


p 


( 21 ) 




The logarithmic derivatives of the structure functions is used to probe intermittenc}^. The logarithmic 


derivatives of Eq. (21), S 4 with respect to 82 ^ for all non-tethered molecules in the turbulent domain, is 
shown in Fig |13a[ The resulting curv e shows similar behaviour to studies of Lagrangian particles in turbulent 
ffow^^. The dip region in Fig |l3a| has been associated with occurrence of small-scale vortex filaments in 
previous turbulence studie^^. The presence of similar behaviour in a molecular simulation may be due to 
rotation inside molecular shells. These rotating vortex like motions can be seen by observing the trajectory 
of a single molecule, as shown in figure [T^ Much of the time is spent orbiting in molecular cages, although 
the meanfiow in x results in a much greater x displacement over time. Due to the low density of the system, 
the molecules will occasionally move large distances. Levy flights, associated with super- diffusion^. This 
behaviour can be seen by looking at the molecular mean squared displacement, defined as. 


1 ^ 


( 22 ) 


i=l 


where the molecules positions are calculated with the x meanfiow removed = ri{t)-\-At [vi{t) — u{t)]. 

The meanfiow, fx, is obtained for a given molecules based on its wall normal position using a spline fit to 
the current velocity profile. The mean square displacement is shown in Fig |14a[ The results from the 
turbulent flow are compared to a periodic equilibrium box and similar simulation of laminar flows with 
matched density p = 0.3 and temperature T « 0.5. Results for a dense fluid, p = 0.81 and T = 0.78, using 
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(a) Logarithmic derivatives of structure function^^, 



(b) Trace of molecule coloured by magnitude of 
particle velocity. 


FIG. 13: Molecular trajectory and structure functions 


the full LJ potential, Vc = 2.25, are also shown. This is consistent with the work of Rahmanl^ who verified 
the simulation diffusion using experimental data. The behaviour is similar for all three systems, although 
despite removing mean velocity, the streamwise diffusion components remain larger. This is attributed to 
the stick-slip behaviour near the walls which results in molecules being dragged at greater velocity than 
predicted from the mean profile. 

In all cases, initially before the molecules collide, they move in the ballistic region with super-diffusive 
behaviour, (r^) oc Eventually, the mean squared displacement becomes linear as predicted by Gaussian 
statistics, (r^) oc t, with the dense fluid case becoming linear earlier than the low density cases. Perhaps 
surprisingly, both the laminar and turbulent flows observe the same mean square displacement as the equi¬ 
librium fluid. It has been shown that diffusion is promoted by strong strain rates in a manner proportional 
to the square root of strain rat^^. As the domain is large and wall speed low, the strain rate is 7 = 0.0035 
and the diffusion enhancement is minimal in both laminar and turbulent cases. The vortices observed in 
the minimal channel seem to have negligible impact on diffusion, which appears to be dominated by the 
molecular structure for low shear rates. 

The autocorrelation of velocity also appears to be dominated by the molecular structure, with a rapid 
decrease seen in the equilibrium, laminar and turbulent flow cases. Fig. |14b| The differences in mean flow 
in the various cases explains the convergence of the autocorrelation to different final values for the various 
cases. The dense fluid decorrelates more rapidly, exhibiting a negative region due to a bounce back after 
collisions. Despite the dominance of molecular scale effects on short term autocorrelation, it appears that 
the hydrodynamics of the vortices may be evident at long times. The insert in Fig. |14b| appears to show 
a further decorrelation in velocity at longer times, possibility due to turbulent mixing Gonflrmation of this 
behaviour would require longer studies to fully explore this phenomenon. However, Lagrangian statistics 
applied to molecular dynamics shows great potential insight into the smallest scale of vortex dynamics in 
turbulent flow. 


IV. DISCUSSION 

In order to focus only on the fluid instability, in this work an incompressible and isothermal GFD solver 
was chosen for comparison. The molecular simulation includes compressibility, temperature and viscosity 
changes as well as thermal fluctuations. A more detailed continuum simulation could certainly include 
these features with compressibility and a changing viscosity coupled to the energy equation. In addition, a 
fluctuating hydrodynamics GFD solver could be employed to reproduce the thermal effects observed in the 
molecular channel. The molecular dynamics minimal channel appears to show remarkably good agreement to 
the continuum model, despite compressibility and temperature change. A lower temperature could be chosen 
or a more aggressive thermostatting mechanism could have been employed to remove this heat. However, the 
increase in temperature due to shear is a unique feature of molecular modeling and it is of interest to include 
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(a) Mean square displacement, insert shows linear (b) Autocorrelation of velocity, insert shows long 

axis near zero and linear (Gaussian) line is shown time for minimal channel flow, 

and labelled. 


FIG. 14: Molecular diffusion and autocorrelation functions. Results from Rahmanl^ (•) with matched LJ 
potential Vc = 2.2, p = 0.81 and T = 0.78 (• • •)? WCA potential at p = 0.3 and T ^ 0.5 for an 
equilibrium case (- -), wall driven Laminar flow (—) and the minimal channel flow (—) 


it in the most realistic manner possible. Despite not being in thermodynamics equilibrium, the breakdown 
and regeneration mechanism is robust enough to repeat twice in the presented work. The turbulent statistics, 
spectral range of energy scales and probability density function are in line with continuum observations. 
The presented study is clearly reproducing many of the features unique to turbulent simulation. By the end 
of the second regeneration cycle, the system energy appears to have reached a stable value, the Reynolds 
number if greater than 400 and the minimal channel flow instability still appears to be present. The rate of 
increase in kinetic energy also appears to be almost zero by the end the second cycle. 

The main feature of this flow, compared to previous NEMD studies, is that it includes time evolving 
hydrodynamics component as well as changing thermodyna mics, resulting in a range of length and time 
scales. The molecular definition of the kinetic pressure tensor, ^{PiPi/rrii)^ and the Reynolds stress tensor, 
pu'u'^ in Eq. (16) are dimensionally and mathematically the same quantity on different timescales. Osborne 
Reynolds actually proposed the decomposition of velocity, u = u^u' ^hy analogy with the kinetic theor}!^. 
Eluctuations which contribute to Reynolds stress simply becomes kinetic pressure once the time and length 
scales of motion become small enough. This is particularly interesting as certain researchers predicts a 
range of scales below the beginning of the dissipation ranged. The range of scales observed in the spectra 
of section III C[ the molecular and averaged probability density function |IIID| and vortex like motions in 
molecular trajectories of section III E| all suggest that molecular studies may provide interesting insight into 
turbulence. 

Reynolds shear stress and shear pressures play similar roles in turbulent and laminar flow respectively. 
Consider an MD simulation of steady state laminar Couette flow. The shear stress from Eq. (14) is given 
by the x momentum flux and force on a plane in the y direction. 


1/ AT AT 

^^xydSy — ~ EE fxij ^^yij 

' i=l i^j 


AT V 

E PxiPyi iq \ _ ^ 

^>^yi 1 — ^laminar 
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(23) 


On average the kinetic pressure and configurational stresses sum to a constant value in laminar flow. The 
equality of Eq. (23) is only true in a time averaged sense, with velocity and stresses fluctuating due to 
what is often termed molecular ‘noise’. In turbulent Couette flow, a similar relation is seen by writing the 
Navier-Stokes equation in terms of Reynolds averaged quantities^. 


YLxydSy pu'v' — Cturbulent 


(24) 


The sum of Reynolds stress and viscous stress is a constant. The form of Eqs. (23) and (24) are strikingly 


similar, both representing a balance between shear effects due to fluctuations and a stress based shearing 
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term. Equation (24) is used to motivate the Boussinesq approximatioi]^ in turbulent modeling which 
introduces an eddy viscosity coefficient to model the details of turbulent fluctuations. Comparison with 
Eq. (23) suggest that molecular detail is essential approximated by the continuum model in the form of a 
viscosity coefficient. Due to a large number of molecules, 0(10^^) molecules per Im^ of air, the molecular 
viscosity is a statistically much better approximation than a turbulence based closure model. 

One possible application of the present work is in the exploration of the transition to turbulence. The effect 
of sub-continuum scales on turbulence is not knowrP^ altho ugh it has been suggested that these microscale 
motions are important to the development of fluid instability^^^. Some authors, even suggest the possibility 
that organized flow may actually originate from the smallest scales, rather than simply dissipate downward 
from the largest scales^. Typically a CED solution requires artificial random noise to trigger the transition 
from laminar to turbulent flovf^. Eor an initially laminar case, the molecular level ‘noise’ ma y p rovide a 
route to transition into a turbulent state as proposed, based on experimental evidence, by Muriell^. Indeed, 
recent work using fluctuating hydrodynamics have shown that thermal fluctuation provide a mechanism 
for the transition to turbulence^. To explore this, a low temperature MD simulation was performed with 
an effective Reynolds number of Re ^ 700 and initial velocity field based on steady state laminar Couette 
solution. The flow remained laminar, suggesting that a larger system (and higher Reynolds number) would 
be required to observe transition from molecular fluctuations in this case. As in a CED simulation at 
Re = 460, the MD solution remains laminar, suggesting that the inherent molecular noi^ at least under 
the current conditions, is not sufficient to induce transition. In recent work, McWhirterl^ suggested that 
the global stability of MD Couette flow should be similar to its continuum counterpart. This could be 
further explored by introducing optimal perturbations to the flow or running the system at higher Reynolds 
numbers. 


Another unique feature of molecular simulation is the ability to explore realistic surface fluid interactions. 
In this study, the wall is solid argon, which forms a perfect face centered cubic crystal structure. Wall-fluid 
interaction, are set to unity, although more realistic walls, along with wall textureJ^ and polymer 
coatings, can be easily constructed in the molecular system. Multiple fluid phases are trivially modelled by 
molecular systems which avoid the discontinuity presented by a three phase contact line and captures bubble 
nucleation. The impact of rapid heating with phase change could be modelled for a turbulent channels, an 
important problem for coolant fluids in micro and nano scale devices. A promising line of research applies 
molecular dynamics for fluid simulation as part of a coupled MD/CED simulation. In coupled simulation, 
an MD so lver is employed in the near wall region and the bulk of the domain is efficiently modelled by CED 
simulatiorP^E^lEIl. As a results, molecular boundaries can be included in higher Reynolds number simulation 
and be used to understand the dynamic interplay of nano-scale boundaries and turbulent flow. This is a 
logical extension of currently employed CED coupling between near wall DNS in LES or RANS solvers to 
provide more detailed modeling at fluid solid interfaces. 

This work has shown that the minimal flow unit, though to be a central component of near wall turbulence, 
can be modelled using a molecular dynamics description. The bursting cycle continues to occur in a similar 
manner to the incompressible continuum model, despite compressibility, significant heating and change of 
viscosity in the molecular model. This suggest that the non-dimensionalised Navier-Stokes equations model 
nonlinear effects which are potential valid to very small system sizes as shown by the good agreement with 
the more fundamental molecular model. The molecular model has great potential to provide insight into 
the nature of turbulent energy cascade below the diffusive range, the fluid-solid interaction and perhaps 
a mechanism for transition from the thermodynamics fluctuations which are an integral part of molecular 
simulation. 


V. CONCLUSION 


A minimal channel planar Couette flow has been simulated using molecular dynamics (MD). The modeled 
MD fluid shows turbulent streaks breakdown and reformation with associated vortex regeneration consistent 
with experiments and computational fluid dynamics (CED) simulations. Reynolds averaged channel statis¬ 
tics show excellent agreement to both an equivalent CED simulation and literature results. A molecular 
form of the law of the wall exposes near wall molecular stacking due to stick slip behavior at the wall. 
The occurrence of the regeneration cycle in MD, which is central to turbulent production, provides strong 
evidence that turbulence like behaviour can be reproduced in a molecular simulation. The insight provided 
by molecular dynamics are explored through velocity spectra, probability density functions and Lagrangian 
statistics. The similarity between the definition of molecular pressure and Reynolds stress suggests a contin- 
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ual cascade where eddying motions don’t simply dissipate to incoherent heat at a certain scale. In this way, 
molecular dynamics is a more fundamental model which requires no assumption of a viscosity coefficient and 
can explicitly models the full energy cascade. As a result, simulation using molecular dynamics has great 
potential in the field of fluid dynamics turbulence research. Simulating turbulence at the molecular scale 
can provide insight into the minimum turbulent eddy and the regeneration mechanism of turbulence from 
molecular origins; explore transition from laminar to turbulent flow; test the limitations of continuum models 
and aid the development of nano-scale modeling methodologies essential in many emerging technologies. 
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